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Doppler spectroscopy has detected 136 planets around nearby stars 1 . A major puz- 
zle is why their orbits are highly eccentric, while all planets in our Solar System are on 
nearly circular orbits, as expected if they formed by accretion processes in a protostellar 
disk. Several mechanisms have been proposed to generate large eccentricities after planet 
formation, but so far there has been little observational evidence to support any particu- 
lar one. Here we report that the current orbital configuration of the three giant planets 
around Upsilon Andromedae 2 3 (v And) provides evidence for a close dynamical interac- 
tion with another planet 4 , now lost from the system. The planets started on nearly circular 
orbits, but chaotic evolution caused the outer planet (v And d) to be perturbed suddenly 
into a higher-eccentricity orbit. The coupled evolution of the system then causes slow 
periodic variations in the eccentricity of the middle planet (v And c). Indeed, we show 
that v And c periodically returns to a very nearly circular state every 9000 years. Our 
analysis shows that strong planet-planet scattering, one of several mechanisms previously 
discussed for increasing orbital eccentricities, must have operated in this system. 

The innermost planet around v And has a very small orbital eccentricity, as expected from 
tidal circularization 5 , but the two outer planets have large eccentricities, both around 0.3. It 
was quickly recognized after their discovery that the gravitational interaction between these 
two planets causes significant eccentricity evolution on secular timescales (~ 10 4 years). In 
particular, in some early solutions, the middle planet appeared to have its eccentricity varying 
periodically with a large amplitude, from a maximum near the present value to a minimum near 
zero 6 . Later fits to the data suggested that the outer two orbits had arguments of pericenter 
nearly equal to within about 10°. This could indicate an "apsidal resonance," in which two 
elliptical orbits oscillate about an aligned configuration with a small libration amplitude 7 . 
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Two possible explanations have been proposed for this peculiar orbital configuration. The 
first is a dynamical mechanism in which a sudden perturbation imparts a finite eccentricity to 
the outer planet 4 . The subsequent secular evolution causes the middle planet's orbit to become 
eccentric and can leave the two orbits either circulating, or librating with a large amplitude. In 
either case, the eccentricity of the middle planet periodically returns to its initial low value. The 
impulsive perturbation of the outer planet would result naturally from planet-planet scattering, 
which was one of the earliest mechanisms proposed for inducing large eccentricities in extraso- 
lar planets 8,9 . In contrast, the second explanation invokes an adiabatic perturbation of the outer 
planet's eccentricity through torques exerted by an exterior gaseous disk 10 . As the eccentricity 
of the middle planet grows on a similarly long timescale, this would leave the system in an 
apsidal resonance by damping the libration amplitude to zero. This is a natural extension of 
more general "migration scenarios" in which the coupling of a planet's orbit to a gaseous disk 
could both increase eccentricities 11 and lead to orbital decay and the formation of planets with 
very short orbital periods 12 . 

The v And system was the first extrasolar multi-planet system ever discovered by Doppler 
spectroscopy. Since the first announcement of the outer two planets in 1999, the California and 
Carnegie Planet Search team has taken over 350 new radial velocity measurements, making v 
And one of the systems with the tightest constraints on orbital parameters (Table 1 ; see also 
supplementary information). The secular evolution of the system is shown in Fig. 1. While the 
eccentricity of d remains always between 0.2 and 0.3, planet c returns periodically to a nearly 
circular orbit with e < 0.01. As shown in Fig. 2, this is now a property of all solutions consistent 
with the data, in contrast to earlier suggestions that it could happen for some solutions 4,6 . 

We have also re-examined the possible presence of apsidal resonance in the system. Re- 
markably, we find that the allowed solutions all lie very close to the boundary between librating 
and circulating configurations (Fig. 3). As a consequence, all librating systems have libration 
amplitudes close to 90°. As shown in Fig. 3, this behavior is confirmed by direct numerical 
integrations of all three planets for a large sample of systems covering the allowed parameter 
space of solutions. 

Our results do not change significantly if the assumption of a coplanar system viewed 
edge-on is relaxed. Inclination angles can be constrained by requiring dynamical stability of 
the system 6,7,16,17 . Using the latest data we find that, for a coplanar configuration of three 
planets with the best-fit orbital parameters, the system becomes dynamically unstable when 
sini < 0.5 (where sini = 1 corresponds to edge-on). If we relax coplanarity, we find that the 
system becomes dynamically unstable whenever the relative inclination is greater than about 
40°. Throughout the full range of allowed inclination angles we find qualitatively similar be- 
havior for the secular evolution of the system. For example, Fig. 2 shows that the eccentricity of 
planet c would still return periodically to a value near zero. The secular evolution of the system 
should be reevaluated if future observations of v And were to discover an additional planet in a 
long period orbit. 

Our analysis clearly confirms that the v And system is evolving exactly as would be ex- 
pected after an impulsive perturbation to v And d 4 . The initial sudden change in eccentricity for 
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the outer planet would be naturally produced by a close encounter with another planet, which 
got ejected from the system as a result. Using our knowledge of planet-planet scattering from 
several previous studies 218-20 , we determined plausible initial conditions for the original, un- 
stable system. The early dynamical evolution of such a system is illustrated in Fig. 4. After a 
brief period of strongly chaotic evolution, lasting ~ 10 3 years, the outer planet is ejected, and 
the remaining two planets are left in a dynamical configuration closely resembling that of v And 
(see supplementary discussion for a more detailed discussion). 

While several other mechanisms (e.g., perturbations in a binary star 22 , resonances 10 , inter- 
action with a gaseous disk 11 ) have been proposed to explain the large eccentricities of extrasolar 
planets, only planet-planet scattering naturally results in an impulsive perturbation, as is nec- 
essary to explain the v And system. All other mechanisms operate on much longer timescales 
and would also affect the eccentricity of planet c, erasing the memory of its initial circular orbit 
(see supplementary discussion for a more detailed discussion). 

Our results have other implications for planet formation. Given the difficulty of forming 
giant planets at small orbital distances, it is generally assumed that the v And planets migrated 
inward to their current locations via interactions with the protoplanetary disk 12 . If this is correct, 
then the small minimum eccentricity of v And c also provides evidence that its eccentricity at 
the end of migration had not grown significantly, in contrast to what some theories predict 11 . 
However, the possibility of formation in situ 23 cannot be excluded by our results. 
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Masses and orbital parameters for the three planets in v And. 



Planet 


P(d) 


e 


u (deg) 


msmi (Mj) 


b 


4.617146(56) 


0.016(11) 




0.6777(79) 


c 


241.32(18) 


0.258(15) 


250.2(4.0) 


1.943(35) 


d 


1301.0(7.0) 


0.279(22) 


287.9(4.8) 


3.943(57) 



Table 1: Results of our new analysis of the v And radial velocity data 2 ' 313 . We have used the 
entire Lick Observatory data set, kindly provided to us by D. Fischer. For conciseness, we 
present only the means and standard deviations on the last two digits (indicated in parenthesis) 
after marginalizing over all other parameters. We list the orbital period (P) in days, the orbital 
eccentricity (e), the argument of pericenter (u) in degrees, and the planet mass times the sine of 
the inclination of the orbital plane to the line of sight (m sini) in units of Jupiter masses (Mj). 
More details on our analysis are presented in the supplementary information. 
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Figure 1 : Secular evolution of the planetary system around v And. The top panel shows the 
semimajor axes (thick lines), as well as the periastron and apastron distances (thin lines), for 
the outer two planets. The lower panel shows the evolution of the orbital eccentricity for each 
planet. Note that both planet c (dashed) and planet d (dotted) have a significant eccentricity at 
the present time (t = 0), but that the eccentricity of c returns periodically to very small values 
near zero. The results shown here were obtained by direct numerical integration using our best- 
fit parameters. All direct iV-body integrations presented in this paper were performed using 
Mercury 14 , version 6.1. 
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Figure 2: Probability distribution for the minimum eccentricity of planet c. We draw initial 
orbital elements for planets b, c, and d from their posterior probability distribution and evolve 
each system according to classical second-order secular perturbation theory. We find that all 
allowed orbital solutions have the eccentricity of planet c oscillating from a maximum value 
slightly larger than the present value to very nearly zero. The solid curve corresponds to copla- 
nar orbits viewed edge-on. The dotted curve shows the result for orbits with random orientations 
to the line of sight, but requiring dynamical stability. This implies relative inclinations < 40°. 
Note that systems with relative inclinations > 140° are also dynamically stable. Although such 
retrograde orbits are unlikely on theoretical grounds, our conclusions are robust to this possibil- 
ity. Since the secular perturbation theory averages over the orbits, it is also valid for retrograde 
orbits. The fact that all allowed solutions result in the eccentricity of planet c returning to a very 
small value can be understood easily from lowest-order secular perturbation theory, 4,15 where 
the eccentricity vector of each planet can be described as the sum of three rotating eigenvectors 
in the (ecoscu, esinu;) plane. The eigenvector representing the effects of planet b on planet c 
has a very small amplitude and can be neglected. For the particular configuration of v And, the 
two dominant eigenvectors describing planet c have very nearly the same length. Depending 
on whether the eigenvector with the faster rotation (higher eigenfrequency) has a slightly larger 
or smaller length than the other, the vector sum will be rotating around 360° (circulation) or 
oscillating with an amplitude close to 90° (libration), respectively. Whenever the two vectors 
are anti-aligned, the magnitude of their sum, i.e., the eccentricity of the planet, is very close to 
zero; when they are aligned, the eccentricity is maximum. 
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Figure 3: Observational constraints on the secular evolution parameters. The eccentricity ratio 
is plotted as a function of the difference between arguments of pericenter for planets c and d, 
all at the present epoch. We show the 1-, 2-, and 3-er contours for the posterior probability dis- 
tribution function marginalized over the remaining fit parameters. The thick contours assume 
that the radial velocity variations are the result of three non-interacting Keplerian orbits viewed 
edge-on, while the dotted contours include the mutual gravitational interactions of the planets 
when fitting to the radial velocity data (only 1- and 2-er contours are shown). The thin solid 
line shows the separatrix between the librating (upper left) and circulating (lower right) solu- 
tions according to the classical second-order perturbation theory (neglecting the inner planet b). 
The dotted lines on either side show the variation in the location of the separatrix due to the 
uncertainty in the remaining orbital elements. The data points show the results of our direct nu- 
merical integrations for the full three-planet system: triangles (squares) indicate that the system 
was found to be librating (circulating). Note that, regardless of the assumptions, the separatrix 
runs right through the \-a contours. We find the median libration amplitude of the librating 
systems to be about 80°. 
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Figure 4: Dynamical evolution of a hypothetical planetary system similar to v And. After a 
brief period of dynamical instability, one planet is ejected, leaving the other two in a config- 
uration very similar to that of v And c and d. The innermost planet was not included, as it 
plays a negligible role. All planets were initially on nearly circular orbits. The initial periods 
of the two inner planets were 241.5 days (equal to the period of v And c) and 2100 days. An 
additional planet (solid), of mass 1.9 Mj, was placed on an orbit of period 3333 days, so that 
the outer two planets were close to their dynamical stability limit 21 . The fourth planet (solid) 
interacted strongly with the third (dotted), while the second planet (dashed) maintained a small 
eccentricity. After the ejection, the two remaining planets evolve secularly just like v And c 
and d (cf. Fig. 1). While this simulation illustrates a case in which v And had only one extra 
planet, our results do not preclude the existence of even more planets at larger distances. In fact, 
the presence of another planet in an even longer-period orbit could be responsible for triggering 
the instability after a long period of seemingly "stable" evolution and help raise more quickly 
the pericenter of the planet which perturbed v And d. Note that the timescale to completely 
eject the outer planet from the system (after ~ 9000 years in this particular simulation) is much 
longer than the timescale of the initial strong scattering. After this initial brief phase of strong 
interaction, the perturbations on the outer planet are too weak to affect significantly the coupled 
secular evolution of v And c and d. Thus, the "initial" eccentricity of v And c for the secular 
evolution is determined by its value at the end of the strong interaction phase, rather that at the 
time of the final ejection. 
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Supplementary Methods 



We model the observations as resulting from three planets on independent Keplerian or- 
bits. We use a Bayesian framework 51 , to constrain the orbital elements and masses with the 
radial velocity observations, as well as the "stellar jitter" (radial velocity variations due to stel- 
lar phenomena such as convection, star spots, and rotation). We assume priors uniform in 
the logarithms of orbital periods, velocity semi-amplitudes, and the stellar jitter. Initially, we 
assume that the orbits are coplanar and viewed edge-on. (Later, we relax this assumption by 
adopting instead an isotropic distribution of inclinations but rejecting configurations that are not 
dynamically stable.) We assume uniform priors for the eccentricities and the remaining angles. 

We use the techniques of Markov chain Monte Carlo (MCMC) and the Metropolis-Hastings 
algorithm with the Gibbs sampler to sample from the posterior probability distribution for the 
masses and orbital parameters' 52 . Since our Bayesian analysis closely follows the methods de- 
veloped for single-planet systems with no stellar jitter by Ford s2 , here we only discuss the 
generalizations to the algorithm which were necessary to apply it to multiple-planet systems 
and to allow for an unknown stellar jitter. Each state of the Markov chain includes the five fit 
parameters for each planet (orbital period, P, velocity semi-amplitude, K, orbital eccentricity, 
e, argument of pericenter, cu, and mean anomaly at the specified epoch, M a ), the unperturbed 
stellar velocity, C, and the magnitude of the stellar jitter, Oj. The jitter is assumed to be Gaus- 
sian and uncorrelated from one observation to the next. We use the Gibbs sampler and Gaussian 
candidate transition functions which are centered on the parameter values from the current state 
in the Markov chain. The scales of the candidate transition functions were chosen based on a 
preliminary Markov chain so as to result in acceptance rates of nearly 40%. The computational 
efficiency of MCMC can also be significantly affected by correlations between variables. To 
improve the computational efficiency we added candidate transition functions which were based 
on several auxiliary variables based on combinations of the fit parameters. These auxiliary vari- 
ables are esincu, ecoscu, uo + M , uo — M Q , P 2 / 3 (l — e), and P 2 / 3 (l + e), for each planet, as well 
as uo b ± uj c , and uj c ± uj d (where cu b , uj c and cod are the arguments of pericenter for planets b, c, 
and d, respectively). The acceptance ratio was determined according to the Metropolis-Hastings 
algorithm to reflect our choice of priors described in the previous paragraph. The use of this 
enlarged set of candidate transition functions significantly improved the mixing and hence effi- 
ciency of our Markov chains. It is important to note that our sampling procedure does not make 
any assumptions about the posterior distribution for the orbital elements or masses. Therefore, 
our approach allows us to take into account correlations between various orbital parameters, in 
contrast to previous simpler analyzes. 

We have computed five Markov chains, each of which contains over one million states. We 
have performed severals checks to verify the reliability of our Markov chains. The acceptance 
rates were between 0.37 and 0.45 for trail states generated by the candidate transition functions 
for each variable or auxiliary variable in each chain. We have also verified that the resulting 
distributions show excellent agreement across all five chains. For example, the Gelman-Rubin 
test statistic, R, approaches 1 from above as a Markov chain approaches convergence. While 
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sets of Markov chains with R less that 1.1, or even 1.2, are frequently used for inference, for the 
Markov chains used in this analysis, R was less than 1.001 for each fit parameter and auxiliary 
variable, and for the predictions of the radial velocity at 40 future epochs. Thus, we are very 
confident in the reliability of our Markov chains. 

We use the results from our Markov chains as the basis for our dynamical analyzes. When 
considering non-coplanar cases, we independently draw the inclinations from isotropic distri- 
butions and the longitude of ascending node from a uniform distribution. For each state in the 
Markov chain, we calculate the planet masses and semi-major axes iteratively and treat them 
as Jacobi elements. Finally, we have performed direct iV-body integrations of systems sampled 
from our Markov chains. For these integrations, we chose a variety of initial epochs varying by 
several times the orbital period of planet d. By comparing the %2 of the fit determined by the 
analytic and iV-body models, we have determined that our results are not affected by including 
the effects of mutual planetary perturbations in the fitting procedure. 



Additional References for Supplementary Methods 
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Supplementary Discussion 



Choice of Initial Conditions for Figure 4 

We have adopted this simple two-planet configuration for computational convenience. In 
this configuration (two planets on circular orbits perturbing each other significantly, with the 
other planets distant enough to be negligible), the dynamical stability limit is known analytically 
and is sharply defined 21 . Therefore we place the two planets very near that limit in our initial 
condition. As discussed in previous papers, a variety of scenarios would lead naturally to two 
planets approaching this limit 81718 . For example, the outer planet might be migrating inward 
through coupling with an outer gaseous disk. Once the stability limit is reached, the system 
evolves quickly (on the orbital timescale) until the strong scattering occurs and one planet is 
ejected, while the gas becomes irrelevant (as the viscous timescale is much longer). 

In addition, the chaotic evolution of a multi-planet system (with more than two planets 
perturbing each other significantly) can also easily lead to strong scattering, sometimes after 
a long period of seemingly "stable" evolution 20 ' 21 . Although equally plausible, this is less 
convenient computationally since, with more than two planets, the stability limit is not known 
analytically and not sharply defined, so numerical experimentation would be needed to find a 
case that could produce a final state resembling the current v And configuration, possibly after 
a very long dynamical integration. 

With more than two planets, the timescale for growth of the instability and the occurrence 
of a strong scattering involving v And d can be arbitrarily long 21,s3 and could easily exceed the 
~ 10 7 yr timescale on which the gas is expected to be lost from the protoplanetary disk. If gas 
were still present when the scattering occurs, the final outcome would be modified, but only 
the details of the dynamical evolution would change. For example, if gas drag produces some 
damping of the eccentricity, then a slightly stronger scattering may be needed to produce the 
same final eccentricity for the retained planet. 

Finally, we want to emphasize that, while reproducing the exact parameters of any particu- 
lar observed system always requires "fine tuning" in an obvious sense, what is important here 
is that our mechanism can — naturally and without any fine tuning — provide the very short 
timescale required for the initial perturbation that left v And d on an eccentric orbit while leav- 
ing v And c on a perfectly circular orbit. 

Perturbations Due to a Gas Disk 

The possibility that the impulsive perturbation to v And d was delivered by a massive exte- 
rior gaseous disk with a large viscosity was mentioned in a previous study 4 . However, we can 
show that this possibility is now firmly excluded for this system. If the eccentricity of v And d 
had been induced by an outer disk, the eccentricity growth time would have to be considerably 
shorter than the secular timescale. In addition, after the eccentricity grew to ~ 0.3, the pertur- 
bation must stop suddenly. Otherwise the eccentricity of the middle planet, v And c, would 
also start growing and its "initial" value would no longer be compatible with the minimum we 
derive in the secular solution. 



13 



Quantitatively, this is a very stringent requirement. We have performed additional numerical 
integrations for the outer two planets, starting with their current masses and orbital periods, but 
on circular orbits. In addition to the mutual gravitational perturbations, we include a simple 
semi-analytic model of angular momentum dissipation acting on v And d only. The dissipation 
rate (J) is constant for a time At and then disappears (completely and instantaneously). We 
have performed multiple simulations holding the product J At fixed to reproduce the current 
eccentricity of v And d. If we impose the constraint that the eccentricity of v And c must 
remain < 0.01 (the current best-fit value of the minimum is 0.005) after At, then this model 
provides an upper limit of At < 100 years. Thus, it is not sufficient to impose an eccentricity 
growth time shorter than the secular period of ~ 10 4 years. 

A timescale for eccentricity growth by viscous coupling to a disk as short as < 100 years 
would require both an implausibly massive disk and a very high effective viscosity. The possi- 
bility of eccentricity growth caused by interactions with a disk is rather controversial 4 , partic- 
ularly for planets less massive than 10 — 20 Mj. Nevertheless, we have estimated the timescale 
for eccentricity using the only detailed theory of eccentricity excitation by viscous coupling 
to a disk in the astrophysical literature 55 . We find that a timescale for eccentricity growth as 
short as ~ 100 years would require a mass in the relevant part of the disk, Sr 2 ~ 40 Mj, ten 
times more than the mass of the planet, even with an implausibly large disk viscosity param- 
eter, a ~ 0.1. Instead, using more typical parameters 55 (for a 1 Mj planet at 1 AU, with disk 
parameters r/h = 25, a = 0.001, and w/r = 0.5), it would require ~ 10 7 years for the eccen- 
tricity to grow from 0.01 to 0.3. Eccentricity growth on this very long timescale would instead 
lead to apsidal resonance between v And c and d with a small amplitude libration 10 , which we 
have ruled out in the main text. In addition, neither this theory 55 nor any other published theory 
of eccentricity excitation due to a gas disk provides a mechanism for making the perturbation 
cease quickly (< 100 years) after a phase of very rapid eccentricity growth. 

What is truly unique about planet-planet scattering is that it provides both a very short 
timescale for eccentricity growth (the dynamical timescale on which the instability develops) 
and the same very short timescale for removing the perturbation (since the extra planet gets 
ejected as a result of the scattering). 
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